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Resonant relaxation and the warp of the stellar disc in the 
Galactic centre 



o 

(N 

> 

o 

a: 

6 

in: 

> 

O 
O 

o 

\6' 
o 
o 



Bence Kocsis^'^'^* and Scott Tremaine^f 

^ Harvard-Smithsonian Center for Astrophysics, 60 Garden St., Cambridge, MA 02138, USA 
^ Institute for Advanced Study, Princeton, NJ 08540, USA 
^ Einstein Fellow 



Received 



ABSTRACT 

Observations of the spatial distribution and kinematics of young stars in the Galactic 
centre can be interpreted as showing that the stars occupy one, or possibly two, discs 
of radii ^ 0.05-0.5 pc. The most prominent ('clockwise') disc exhibits a strong warp: 
the normals to the mean orbital planes in the inner and outer third of the disc differ 
by ~ 60°. Using an analytical model based on Laplace-Lagrange theory, we show that 
such warps arise naturally and inevitably through vector resonant relaxation between 
the disc and the surrounding old stellar cluster. 
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1 INTRODUCTION 



Black holes of mass are found in the centres 

of most galaxies. These exotic objects are the engines that 
drive quasars and other active galactic nuclei, and may play 
an important role in galaxy formation through feedback to 
the interstellar medium. Our own Galaxy contains a central 
black hole of mass ~ 4 x 10® M0 associated with the radio 
source Sgr A*, and thus provides a unique opportunity to 
explore the interactions of a nuclear black hole with the 
surrounding gas and stars, at spatial resolution far greater 
than can be achieved for any other galaxy. 

Among the more remarkable components of the Galac- 
tic nucleus is the population of ~ 100 massiveyoung stars 
found in the central parsec (Ipc = 26arcsecjj. These are 
mostly O supergiants and Wolf-Rayet (WR) stars, with 
masses >2OM0, formed in a burst lasting <2 Myr ap- 
proximately 6 ± 2 Myr ago (|Paumard et al.l l2006l ). Proper 
motions an d radial velocities are available for most of 
these stars (jBartko et al.l [2OO 9I). In the standard descrip- 
tion, about half of the massive stars between 1 arcsec and 
10 arcsec belong to a rotati ng system (the 'clockwise disc', 
iLevin fc Beloborodovl [20031 ) . which can be modelled as a 
warped disc with local thickness (standard deviation of the 
inclinations) of abo ut 14° and mean eccentricity of 0.3-0.45 
ijBartko et al.ll2009l ). The warp is substantial: the symme- 
try axis of the disc varies by 60°-70° between 1 arcsec and 
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^ We exclude the separate population of stars ('S-stars') found 
at much smaller radii, <J0. 5 arcsec ~ 0.02 pc. 



10 arcsec radius. About 20% of the stars appear to be mem- 
bers of a second rotating system (the 'counter-clockwise 
disc'), which is thicker than the clockwise disc and inclined 
by 115° to the clockwise disc in the same radius range. The 
properties and even the existence o f the counter-clockwise 
disc are controversial (|Lu et al.ll2009|). even though its s tatis- 
tical significance is 7a according to Bartko et al. (I2OI0I). The 
total mass of the two discs is 5-10 x 10^ M© ( Paumard et al.l 
|2006|). Inside 1 arcsec — 0.04 pc there is a sharp cutoff in the 
density of WR/0 stars. The surface density of the clockwise 
disc is 



E cx r 



(1) 



with S :^ 1.4 ± 0.2 between 1 arcsec and 15 arcsec 
jBartko et al.ll2010l ). T he disc(s) are surroun ded by a spher- 
ical cluster of old stars jBuchholz et al.ll2009l ). The cluster of 
old stars is much more massive than the disc(s) - 5 x 10^ M0 
inside the outer edge of the disc at 10 arcsec - but still <10% 
of the mass of the central black hole, so the disc(s) are nearly 
Keplerian. 

The disc(s) present a number of puzzles: 

• The existence of young, massive stars implies that star 
formation must have occurred recently in the central parsec. 
This is surprising, since tidal forces from the black hole sup- 
press star formation unless the density of the protostellar 
cloud is orders of magnitude larger than cur rently observed 
in this region (|Morridll993l : lAlexandejIioOsI ) . 

• The presence of two distinct discs suggests that there 
were two separate star-formation events. But then why do 
the stars in t he two discs appear t o have the same age to 
within IMyr (|Paumard et al.ll2006l )? 

• The complexity of the dynamical models (two intersect- 
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ing discs, warps, etc.) needed to explain the kinematic and 
spatial distribution of the disc stars suggests that some other 
structure may provide a better description of the data. What 
is the nature of this structure and why is the distribution of 
young stars so complicated? 

In this paper we shall focus on one aspect of the 
Galactic-centre disc(s), their substantial warp. We shall ar- 
gue that the warp arises from slowly varying stochastic 
torques exerted on the disc by the surrounding cluster of 
old stars, through the process sometimes called vector res- 
onant relaxation. Other properties of the disc(s) may also 
arise through resonant relaxation, a discussion that we de- 
fer to future papers. 

There is an extensive literature on warps of stel- 
lar and gasous discs i n the galacti c context (for gen- 
eral reviews see, e.g., iBinnevI [l99^: iNelson fc Tremaind 



1 19961 : iBinnev fc Tremaind 20081 ) . In subparsec scale ac- 



cretion discs, warps can arise through relativistic frame 
dragging (L ense-Thirring precess i on) a round a spinning 
black hole ([Bardeeii fc PettersonI 1 19751 ). radiation pres- 
sure (iPettersonT 19771 : Pringld 19961 ) . gravitational torques 
due to massive tori such as the molecular torus i n 
the Galactic centre (|Navakshinl l2005l : ISubr et al. 

I l2009l ). 

a binary companion s uch as a second black hole or- 
biting inside the disc llPap aloizou. Terguem. fc LinI 1 19981 : 
lYu fc Tremain3l2003l : IYu et al. 2007 ), or stochastic torques 
from a surrounding star cluster l|Bregman fc Alexander! 
I2OO9I ). The self- gravity of a stellar disc can play 
an im portant role in determining the shape of the 
warp faunter fc T oonir? ' 1 9691 : iNelson fc Tremaind 1 19961 : 
lUlubav-Siddiki et al. 2009). Using A'^-body simulations of 
an isolated, initially flat and thin stellar disc t h at re - 
sembles the Galactic-centre disc, ICuadra et al.l (|2008l ) 
showed that the observed warp c annot arise from in- 
teract ions among the disc starj ^ . iNava kshin fc Cuadral 
l|2005l ) and lHobbs fc Navakshiril (|2CI09l ~ have suggested that 
the Galactic-centre disc(s) could have been formed in a 
high-inclination collision betwe en two massive gas clouds. 
iLockmann fc Baumgard^ (|2009l ) have examined the inter- 
action of the stellar disc with a possible massive inclined 
second stellar disc, and showed that the discs dissolve due 
to the mutual torques on timescales comparable to their age. 

In this paper, we examine the evolution of an initially 
thin, flat disc in a near-Keplerian gravitational potential, 
accounting for both the self-gravity of the disc and stochas- 
tic gravitational torques from a surrounding cluster of stars. 
We argue that the most important torques are those that 
arise from the mass distribution of the cluster stars after 
averaging over orbital phase and apsidal angle (vector res- 
onant relaxation) and that the self-gravity of the disc sup- 
presses the excitation of small-scale normal modes so that 
vector resonant relaxation warps the disc, rather than thick- 
ening it. In ^ we discuss the timescales of various pro- 
cesses in the Galactic centre, and establish that vector res- 
onant relaxation with the old cluster stars is significant for 
the Galactic-centre disc(s), whereas most other dynamical 



^ This result is consistent with the analytic arguments below that 
two-body and resonant relaxation among the disc stars is unim- 
portant. 



relaxation processes (e.g., scalar resonant relajcation, two- 
body relaxation, etc.) are not. In S|3] we derive an analytic 
solution to the time evolution of an initially flat stellar disc, 
based on the following approximations: (i) the orbital period 
and the apsidal precession period are much shorter than the 
timescale for the warping of the disc, so we can average 
the motion over the orbital phase and apsidal angle; (ii) ex- 
ternal torques on these timescales are generated exclusively 
by the orbit-averaged mass distribution of the cluster stars 
(i.e., vector resonant relaxation), (iii) the eccentricities of 
the disc stars are negligible; (iv) the warping angle is small 
(in particular, the relative inclination between any two disc 
stars is small compared to their fractional difference in semi- 
major axis); (v) the orbits of cluster stars are uncorrelated 
and independent of the disc (i.e., the back- reaction of the 
disc on the cluster is negligible); (vi) the cluster is spherical 
on average (i.e., the non-spherical component of the clus- 
ter mass distribution is due solely to Poisson fluctuations 
from individual stars). In this case, the problem reduces to 
Laplace-Lagrange theory, in which the secular evolution of 
the disc is described by a quadratic Hamiltonian, and the 
disc is equivalent to a system of point masses interconnected 
with springs and driven by the external forces from the clus- 
ter. This system is integrable as the disc can be decomposed 
into independent normal modes (i.e., independent harmonic 
oscillators). In 521 we derive the stochastic evolution of the 
normal mode amplitudes and consider applications to the 
discs in the Galactic centre. Our conclusions are discussed in 
[J6l The statistical equilibrium of an isolated self-gravitating 
stellar disc is presented in the Appendix. 

In future work we shall present more general numeri- 
cal models that do not require the approximation that the 
inclinations and eccentricities are small. 



2 TIME-SCALES IN THE GALACTIC CENTRE 

We now ask which dynamical processes can shape the prop- 
erties of the disc over its lifetime. 

Two recent estimates of the distance and mass of the 
black hole in the Galactic ce ntre are Rp = 8. ± 0.6 kpc, 
M. = (4.1 ± 0.6) X 10^ Mq (|Ghez et al.l l2008h an d Rp = 
8.3 ± 0.4 kpc, M. = (4.3 ± 0.4) x IQ^Mq (Gill essen et all 
l2009t ). For simplicity we shall adopt _Ro = 8 kpc and M, = 
4x 10^ Mq. At this distance 1 pc = 25.8 arcsec and 1 arcsec = 
0.0388 pc. 

The times and distances derived below are plotted in 
Figure [TJ which also shows the disc age (6 ± 2 Myr) and 
radial extent (0.04-0.6 pc) as a hatched bar. 

Gravitational radius The gravitational radius of the 
black hole is 



^. = ^^^1.18xl0i^cm. 



(2) 



marked in the figure by a dashed vertical line. 



Orbital frequency The orbital frequency Q, of a. star with 
semi-major axis a is given by 



\GM. 



1/2 



= 236 yr 



0.1 pc 



3/2 



(3) 
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Figure 1. Time-scales in the central parsec of the Galaxy. The dashed black lines show the gravitational radius of the central black hole 
(eq. [2J and the age of the Galaxy. The hatched black line shows the age and radial extent of the discs(s) of WR/O stars. The slanted 
black line at the lower right shows the inverse orbital frequency (eq. [Sjl. The magenta lines show the apsidal precession times due 

to general relativity IjJq^ (eq. |4]| and the stellar cluster |a;jv| — ^ (eg. I12II . and the magenta points show the combined apsidal precession 
time |a;GR-|-£jjv|~^. The solid cyan line slanting up to the right shows the Lense-Thirring precession time (eq.|6ll. The red line shows the 
two-body relaxation time for the stellar cluster assuming that m2 = {vn? ) / (m) = IMq (eq. I13|l , and a magenta line shows the two-body 
relaxation time within the disc, assuming m2 = 10 Mq and a disc mass of 5000 Mq (ea. l20ll . The dotted lines show the scalar and vector 
resonant relaxation time-scales in cyan and green, for m2 = 10 Mq (eqs. I24l and l26ll . The solid blue lines show the collision time in the 
stellar cluster for stars with mass m = M0 and radius = Rq and in the disc for m = 20 Mq and radius r* = IORq fcas. l2ll and 12311 . 
Finally, the slanted dashed black line shows the precession time due to the molecular torus. 



The characteristic orbital time Q, ^ (orbital period divided 
by 2tt) is plotted in Figure [T] 

Relativistic precession The apsidal precession rate due 
to general relativity cjgr is given by 



(1 



3(GAf.)»/2 



= 4.11 X 10^yr(l - e^) 



0.1 pc 



5/2 



(4) 

where e is the eccentricity. The characteristic precession time 
oJq^ for nearly circular orbits (e ~ 0) is plotted in magenta 
in Figure [1] Over the radius range of the disc(s), the rel- 
ativistic precession time (|4]) is much larger than the New- 
tonian precession time (|12|) . also plotted in magenta. Thus 
relativistic apsidal precession is unimportant for the disc(s), 
although it is likely to dominate over Newtonian precession 
inside 0.1 arcsec. 



The orbit-averaged Lense -Thirring precession of a star 
with angular momentum L is iLandau fc Lifshit dlioOTl ') 



dL 

It 



_2G2Mts_ 

a3c3(l-e2)3/2 



n,xL = ujuT n, xL, 



(5) 



where n, is the unit vector parallel to the spin axis of the 
black hole and ^ s < 1 is the spin parameter, that is, the 
spin angular momentum of the black hole is sGM"^ /c. The 
characteristic precession time 



4.45 X 10^° yr ,^ 2 3/2 
= (1 - e ) ' 



0.1 pc 



(6) 



is plotted in Figure[T]in cyan, for circular orbits and a maxi- 
mally spinning black hole (s = 1). Lense-Thirring precession 
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is negligible for the disc(s), at least at the radii where they 
are currently observed (^larcsec). 

The star cluster The black hole is surrounded by an ap- 
proximately spherical cluster of old stars. The distribution 
of mass in the cluster can be measured from star counts by 
assuming that the mass density is proportional to the num- 
ber density. The proportionality constant can be estimated 
from the kinematics (radial velocities and proper motions) 
of stars at radii ~ 1 pc, where the stellar mass begins to 
make a substantial contribution to the overall gravitational 
force field. 

Using data from ISchodel et all (|2007l ') , iLockmann et al.l 

l|2009l ) estimate that the mass density at radius r is given 
by 



UN 



p(r) = (2.8 ± 1.3) X 10*^ Mq pc" 



0.22 pc 



(7) 



with 7 = 1.2 inside 0.22 pc and 7 = 1.75 outside 0.22 pc. 
The corresponding enclosed mass is 



0.50 X 10^(r/0.1pc)^■^ 
1.12 X 10^(r/0.1pc)i-2^ 



(8) 



0.92 X 10" 



r < 0.22 pc, 
r > 0.22 pc, 



with an uncertainty of about 50% . 

For comparison, iTrippe et al.l (|2008l ) give 



p(r) = 2.1 X IO^Mqpc"^ 



which yields an enclosed mass 



IH- (r/0.34pc)2' 



i^ = 1.0xlO«Me/ 



0.34 pc 



(9) 



(10) 



where f{x) = x —- tan~^a::. The two expressions for M(r) 
agree to within a factor of two for r^O.3 pc but differ by more 
than an order of magnitude for r <0.05 pc. At r = 1 pc the 
two expressions give M(r) = 1.9 x 10® and 1.7 x 10*^ M© 
respectively, somewhat larger than the independent estimate 
of 1.1-1.5 X lO'' M0 given bv lSchodel eFail l|2009l ). For our 
numerical results we shall use the parametrization ((8|. 

Velocity dispersion Assuming that the velocity disper- 
sion tensor of the stellar cluster is approximately isotropic, 
and solving the Jeans equation for the one-dimensional 
veloc ity dispersion a"(r) (eq. 4.216 in iBinnev fc Tremaind 
l2008l ). we find 



r(r) = < 



280 km s^^ \/0.1pc/r \/l- 0.035 (r/O.lpc)^ ^ 

if r < 0.22 pc, 

if r > 0.22 pc. 

(11) 



250km s-VO-lpc/r, 



Newtonian precession The apsidal precession rate wn 
due to the gravitational field from a spherical star cluster 
is always negative, that is, the line of apsides pr ecesses in 
the o pposite direction to the orbital motion (e.g.. FlYemainel 
I2OO5I ). For a cluster with density p{r) and mass Mir) <C M, , 
the precession rate of a nearly circular orbit with radius r is 



-2TiGp{r)/n{r), oiEl 

2.1 X 10''yr(r/0.1pc)-''^ 
1.3 X 10■*yr(r/0.1pc)''■2^ 



if r < 0.22 pc 
if r > 0.22 pc. 



(12) 

plotted as a magenta line in Figure [T] We also plot as ma- 
genta dots the total precession rate cj = cjgr + ^n- Over 
the radial extent of the disc(s), the characteristic precession 
time is dominated by Newtonian effects and equal to a 
few times 10* yr, more than two orders of magnitude shorter 
than the disc age. 

Two-body (non-resonant) relaxation The two-body 
relaxation tim e for the cluster of ol d ^ars is given by equa- 
tion (7.106) in lBinnev fc Tremaind (|2008i '). 

trelax = 0.34 



G^pm2 In A 



= < 



' 3.6 X loVr (r/O.lpc)"°-^(15/lnA)(M0/m2) 

X fl - 0.035 (r/0.1 pcf-^] ^'"^ , if r < 0.22 pc 

1.7 X lO^r (r/0.1 pc)''-^^ (15/ In A)( Mo/ma), 

ifr>0.22pc, 
(13) 

where InA ~ ln(Af,/m) ~ 15, and the effective mass m2 = 
(m^) /(m) is the ratio of the mean-square stellar mass to the 
mean stellar mass; we shall call this the effective mass. The 
relaxation time treiax for In A = 15 and m2 = is plotted 
in red in Figure [1] 

Unfortunately the effective mass is quite uncertain. We 
assume a stellar mass function of the form 



dnoc ra °'dm for rrin 

then if mmin <^ mmax, 

r (a-2)(Q-3)-Vn.ir 

m2=l (a-2)(3^Q)-' 
[ (2-a)(3-a)-^ 



(14) 



3-a a-2 
'"-max "^min 



a > 3; 

2 < Q < 3; (15) 
a < 2. 



Thus for a standard Salpeter mass function (a = 2.35) with 
ffl-m ax = 100 Ma and m mm = 0.1 M©, m2 = 4.8 M©. For 
the lKroupa et al.l l| 19931 ) solar neighborhood mass function, 
7712 = 0.66 M©. However, there is evidence that the mass 
function in the Galactic ce ntre is much more to p-heavy than 
in the solar neighborhood. iBartko et all l|2010l ) suggest that 
the Galactic-centre disc(s) have a ~ 0.45, in which case 
7712 = 0.677iniax ~ 60 M© . Eveu if the initial mass func- 
tion were known, there are other complications. First, the 
initial-final mass function - the relation between the main- 
sequence stellar mass and the mass of the compact object 
that remains after stellar evolution is complete - is poorly 
known. Second, star clusters, massive gas clouds, and other 
density inhomogeneities also contribute to m2 and may lead 
to values of m2 that ar e orders of mag nitude larger than 
the typical stellar mass (IPerets et al.l [20071. An example in 
the Galactic Centre is IRS 13E, a compact cluster of three 
bright blue supergiants and many fainter components with a 



3 Tills formula was derived for eccentricity e <g 1, but works 
fairly well for e^l. For example, if the cluster density is a power 
law in radius, p{r) oc with 1 < 7 < 3, the error is less than 
30% for all eccentricities e < 0.7. 
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velocity dispersio n suggesting a central mass of a few 10* Mq 
l|Fritz et al.ll2010l ). 

A further uncertainty is that two-body relaxation may 
lead to mass segregation so that the effective mass depends 
on radius. In particular, the relaxation time inside a few 
tenths of a parse c may be dominated by stellar remnants 
rather than stars (jO'Learv et al.|[2009l l. If the mass function 
is broad and dominated by light stars, then the heavy stars 
develop a much steeper cusp with p oc to . This 
affects m2 primarily very close to the black hole, increasing 
1 712 by a factor of ~ 2 at 0.4 pc and a factor ~ 4 at 0.04 pc 
j Alexander fc Hopmaiill2009l : iKeshet et al.ll2009l '). 

The two-body relax;ation time for m2 = 1 M© shown in 
Figure [1] is never shorter than a few Gyr so even if m2 ~ 
10^ M0 two-body relaxation is negligible over the age of the 
disc stars, ^10 Myr. However, two-body relaxation may have 
a substantial effect on old stars in the surrounding cluster. 

For massive stars such as the WR/O stars observed in 
the disc(s), dynamical friction from the star cluster can act 
on a shorter time-scale than troiax- The time-scale for orbital 
decay of a star of mass m on a near-Keplerian circular orbit 
of radius r is given by 

ti,l = ~^S^n^^\nKg{X), (16) 
^^'^ r dt M. M, ^ ^ 

where g{X) = erf(X) - Xerf '(X) and X ^ v/{V2a) where 
V is the circular speed at r l|Binnev fc Tremaine .20081 ^. If 
the density varies as a power law with radius, p oc r then 
X = [^(1 + 7)]^/^. Using u = (GM./r)^^^ and the dispersion 
profile pi|l . we find that g{X) varies from 0.47 for r <^ 0.2 pc 
to 0.57 for r ^ 0.2 pc so for our purposes it is adequate to 
use a constant value g{X) — 0.5. With this approximation 

1.7 X 10** yr (r/0.1 pc)""'^ (12/ In A)(20 Mo/m), 

if r < 0.22 pc, 

1.1 X 10* yr (r/0.1 pc)"-^^ (12/ lnA)(2OM0/m), 

if r > 0.22 pc. 



^fric 



(17) 

We conclude that the action of dynamical friction from the 
star cluster on the disc stars is negligible over the disc age 
of 6 Myr, even for stars as massive as 100 M0 . 

The effects of two-body relaxation between the disc 
stars are uncertain because the properties of the disc(s) are 
poorly determined and the relaxation time is a strong func- 
tion of the root mean squared (rms) eccentricity and inclina- 
tion. For our purposes it is sufficient to approximate the disc 
as a single population of stars of mass m, a nd in this case 
the e ccentricity relaxation time is given by jStewart fc Idal 
I2OOOI ) 



1 d{e^) 
{e2) dt 



-4.5 



In A. 



(18) 



Here (e^) is the mean-square eccentricity of the disc stars, 
E(r) is the surface density, and A ~ {e^)^^^M,/m. This for- 
mula assumes that the rms inclination {i'^)^^'^ ~ 0.5{e^)^''^, 
a typical value seen in relaxed discs and consistent with ob- 
servations of the clockwise disc; and that A ^ 1 (the disc is 



* This result assumes that the distribution function is Maxwellian 
but should be approximately valid for more realistic distribution 
functions. 



'dispersion-dominated'). The inclination relaxes at a rate 
J_d{i^^J^_ (19) 

(z*^) dt ^rclax 

We assume that the surface-density distribution in the 
disc is given by equation ^ and parametrize the disc by 
its total mass Mdisc , which is probably about 5000 Mq 
l|Paumard et al.ll2006l ). Then 



= 4 X 10° yr 



8 (e^)^ 5000 Mq 10 M0 /O.lpc 



(0.3)4 Mii, 



m2 



9 



In A' 
(20) 

where InA = 9.3 corresponds to m = 10 Mq, (e^)^^^ — 0.3. 
Given these parameters, the disc relaxation time is plotted 
in magenta in Figure[T] Even for a large effective mass, m2 ~ 
10^ M0, the minimum relaxation time over the radial range 
of the disc(s) is ^3 x 10^ yr, so two-body relaxation should be 
negligible for most stars over the disc age of 6 Myr. However, 
if the eccentricities and inclinations in the disc have been 
over-estimated, perhaps because of unrecognized systematic 
errors, the scaling troiax ~ e"* implies that the relaxation 
time could be much shorter. 

Physical collisions For a population of stars of a single 
mass, the rate of physical collisions is (|Binnev fc Tremaind 
I200a eq. 7.195) 



coll 



16v^narJ(l -f- 9) 



(21) 



where n — p/m is the number density of stars and B = 
^vl/a^ where = ^2Gm/r* is the escape speed from the 
surface of the star, radius r*. The collision time of stars in 
the cluster is shown in blue in Figure [TJ assuming that the 
cluster is mainly composed of solar-type stars with m = M0 
and Vi, = 618 km s~^. At radii ^0.01 pc, where the collision 
time is less than a few Gyr, collisions are likely to play a 
dominant role in determining the distribution of old solar- 
type stars and the giants into which they eventually evolve, 
which have much shorter lifetimes but much larger radii. 
CoUisional destruction of red-giant envelopes may be respon- 
sible for the depletion in luminous red giants observed inside 
~ 1 pc and the disappearance of the CO spectral feature as- 
sociated with red gian ts in the integrated light (|Alexandeij 
l2005l : lDale et al.]|2009l ). However, collisions with stars in the 
old cluster do not play a major role in the evolution of the 
early- type stars in the disc(s), since the collision time is 
much longer than the disc age, even after accounting for 
the larger radii and masses of the WR/O stars. 

In a dispersion-dominated disc composed of identical 
stars of mass m and r adius r*, the collision rate is (e.g., 
iHeng fc Tremaind I2OO9I ) 

U = 16^f^r.^(0.69 + 1.526), e = ^^^. (22) 



With the same parameters used to derive equation (|20|) . we 
find 



4.1 X 10" yr 5000 Mc 



1 + 2.209 Mdis 



9 = 24.6- 



20 M<7 



^■VlO7?0^ ^ 



0.1 pcy V 
(0.3)^ 



C23~) 

■2OM0 O.lpc {e2) ■ ^ ' 

This result is shown as a blue line in Figure [T] The collision 
time is longer than the disc age. 
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Scalar resonant r elaxation Resonant relaxation 
l|Rauch fc Tremaine|[l996l ) arises from the forces due to the 
orbit-averaged mass distribution of the stars. Because these 
forces vary only slowly, they affect the angular momen- 
tum but not the energy of stellar orbits. Scalar resonant 
relaxation affects the magnitude of the angular momentum 
and thus the eccentricity, while vector resonant relaxation 
affects the direction of the angular momentum vector but 
not its magnitude. The sc alar resonant relaxation time in a 
spherical stellar system is l|Hopman fc Alexanderlliood ) 



Adopting /3„ = 1.83 ±0.03 from lEilon et all (|2009l '). we have 



I31Q? M{r)m2 ' 



(24) 



here u is the apsidal precession rate, the sum of the (neg- 
ative) Newtonian rate ojn (eq. I12p and the (positive) rela- 
tivistic rate ojqr (eq. H]). The di mensionlc ss coefficient Ps 
is estimated to be 1.05 ± 0.02 bv lEilon et~al.l ((20091). Using 
these parameters and effective mass m2 = 10 M0 we plot 
the scalar resonant relaxation time using cyan dots in Fig- 
ure [T] The cuspy minimum near r — 0.007 pc arises where 
relativistic precession cancels Newtonian precession so the 
total precession rate vanished 

The scalar resonant relaxation time is less than 10^" yr 
throughout the central parsec. Thus the eccentricity distri- 
bution of old stars in this region should be relaxed. However, 
the scalar resonant relaxation time exceeds 10* yr through- 
out the radial extent of the disc(s), so eccentricity relax- 
ation of the disc stars is likely to be small. However, it is 
likely that the eccentricity relaxes faster for the most mas- 
sive stars in the cluster, similar to how dynamical friction 
( 161) is faster than two-b ody relaxation for massive stars 



Ranch fc Tremain3 1 19961 ). Assuming that the analogous 



resonant dynamical friction time-scale is inversely propor- 
tional to stellar masfl it is still comparable to the disc age 
only near its inner edge. 

Scalar resonant relaxation can also occur among the 
stars in the di sc(s). The rate f or this process is more difficult 
to determine (|Tremainelll998h . 



Vector resonant relaxation This process affects the ori- 
entation of the angular-momentum vector but not its mag- 
nitude. In a stellar system in which the time-averaged grav- 
itational field is sp herically symmet ric, the vector resonant 
relaxation time is (|Eilon et aklbOO^ fl 



27r 1 



(25) 



^ The sharpness of this minimum is artificial, since the total pre- 
cession rate vanishes at different semi-major ax es for orbits of 
different eccentricities llGiirkan fc Hopmanl^2007^ . 
^ The rate of resonant dynamical friction can be computed us- 
ing the formalism for ordin ary dynamical friction derived by 
iTremaine fc Weinberg l ll984) . in which the friction arises from 
stars that are in near-resonance with an orbiting massive body; 
resonant friction is stronger than ordinary friction because some 
of the resonant denominators are near zero. 

We have set Eilon et al.'s parameter to unity. 



= < 



7.9 X 10® yr ^ Mq /m2 (r /0. 1 pc)°■^ 

if r < 0.22 pc, 
5.3 X 10® yr ^M0/m2(r/O.l pc)°-^ 
X [1 - 0.82(r/0.1 pc)-i-25] if r > 0.22 pc. 

(26) 

This result is plotted in green in Figure [T] for 7712 — 10 M0 . 
At all radii, vector resonant relaxation is substantially faster 
than scalar resonant relaxation. The inclination distribution 
of old stars should be relaxed throughout the stellar cluster. 

Equation (|26p is only valid if the nodal precession rate is 
dominated by the stochastic component of the non-spherical 
gravitational force. Hence this equation may overestimate 
the relaxation rate if the stellar system is not spherically 
symmetric. For example, the non-spherical field from the 
disc is smaller than the stochastic field from the cluster stars 
only if Mdisc<[M(r)m2]^''^ For M(r) ~ 1 x lO'^ Mq (cf. eq. 
[8| and m2 ~ 100 Mq (see discussion following eg. I15p . this 
requires M^isc^S x 10^ Mq, compared to an estimated mass 
of 5 X 10'^ Mq. Thus the approximation that the precession 
rate is dominated by stochastic forces is suspect, and should 
be improved in future models. 

Vector resonant relaxation may be responsible for warps 
in gaseous accretion discs in the centres of galaxies, in partic- 
ular f or the warp of the maser disc in the galaxy NGC 4258 
(Brcg man fc Alexander! l200d V It may also be the principal 
mechanism that i sotropizes the inclinations of S-stars in the 
Gala ctic centre (|Hopman fc Alexander! 12006! : iPerets et al! 
l2009h . In this paper, we investigate whether vector resonant 
relaxation leads to the warped structure of the stellar disc 
in the Galactic centre (|Bartko et al.![200^ ). 

Vector resonant relaxation can also occur among the 
stars in the disc. Nodal precession within the disc is deter- 
mined by the mean gravitational field of the flattened mass 
distribution, rather than the random torques from individ- 
ual stars. The typical nodal precession rate for a star in the 
disc iflfl 

n Mdisc 



(27) 



We then estimate the vector resonant relaxation time from 
equation (I24|) by replacing the apsidal precession rate u with 
the nodal precession rate 1/ and the cluster mass Af (r) with 
the disc mass Mdisc: 

M. 



^rr ,v 



47r 



i2.4 X 10^ yr- 



14° 20 Mp 



i2)l/2 



m2 



0.1 pc 



3/2 



(28) 



This is much longer than both the age of the disc(s) and 
the vector resonant relaxation time-scale due to the clus- 
ter stars. Thus vector resonant relaxation among the disc 
stars is unlikely to play an important role in determining 
the properties of the Galactic-centre disc(sjf|. 



° Note that the rate diverges as the thickness of the disc ap- 
proaches zero, because the torque exerted on a ring by a massive, 
razor-thin disc approaches infinity near the disc. 
^ Nevertheless, it is worthwhile to understand the effects of vector 
resonant relaxation on an old, isolated disc, and this is the subject 
of the Appendix. 
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Warping by a massive perturber The disc(s) can 
be warped by a massive perturber, su ch as the molecu- 
lar t orus outside the disc at 1.5-7 pc (jChristopher et aP 
or an intermediate-mass black hole inside the disc 
l Yu et al] |2007| ). This process was first investigated by 



Laplace in the context of planetary satel lites and has been 
studied by many authors since then (iHunter &: Toomrd 
11969 1: 'Nelson & Trcmainc 'l996'; 'Ulubay-Siddiki ct al."2009); 
seelNay akshi n (2005 ), Lockmann fc Baumg ardt (2009), and 
ISubr et alj ( '2OO9I ) in the context of the Galactic-centre 
disc(s). The precession rate of a circular ring at radius r 
due to a circular ring of mass nip at radius with relative 
inclination I is given by equation (|73|l : i n the limit wher e 
r <C or r ^ Tp the result simplifies to (|Navakshinll2005l ) 



4 M. 



rr 



cos 7, 



(29) 



where r< — min(r, rp), r> = max(r, rp). This formula as- 
sumes that the mass of the perturbed ring is sufficiently 
small that its angular momentum is much less than the an- 
gular momentum of the perturber. 

The characteristic precession time tp = l/vp for an ex- 
ternal perturber (r <^ rp) is then 



ip 



8.5 X 



lOVr^O' 



Mr. 



1.5pc 



0.1 pc 



3/2 



0.5 



cos / 



(30) 

The reference masse s and radii have been s et to resemble 
the molecular torus l|Christopher et al]|2005l ). For these pa- 
rameters, the precession time-scale is longer than the age of 
the disc at r;^0.1pc. However, for r^O.lpc, the molecular 
torus can warp the Galactic-centre disc(s) significantly, at 
least if cos / is not too small (Fig.[T]). Curiously, the molecu- 
lar torus ts nearly orthogonal to the mean orientation of the 
clockwise stellar disc (cos / ~ 0.00 ±0.03). More precisely it 
is nearly orthogonal to the inner parts (cos / ~ 0.05±0.03 for 
r;^0.15pc), but not the outer parts (cos / ~ 0.52 and 0.24 
between r = 0.1 3-0.27 pc and 0.27-0.46 pc, respectively). 
ISubr et al. I (|2009t ) argued that the molecular torus could lead 
to the thickening and warping of the stellar disc if it were 
initially not orthogonal. 

The precession time for an internal perturber (r S> rp) 

is 



5.7 X 



^ . 5000 Mq 
10 yr — 

Mp 



0.1 pc 



0.2 pc 



7/2 



0.5 



cos / 
(31) 

The reference masses and r adii have been set to resemble 
the counter-clockwise disc (|Bartko et al.l [20091 ). For these 
parameters the precessi on time-scale is c omparable to the 
age of the disc; in fact, iNavakshin et al.] (|2006 ) have used 
Ai'-body simulations to set an upper limit of 5 x 10^ Mq on 
the mass of the counter-clockwise disc from the requirement 
that it does not warp the clockwise disc too much. 

We have seen that a rich set of internal and external dy- 
namical processes can affect the evolution and current state 
of the Galactic centre disc(s). A useful first approximation is 
to neglect external tidal fields, two-body relaxation, dynam- 
ical friction, scalar resonant relaxation, and scalar resonant 
friction with the surrounding old stellar cluster, as well as 
two-body and resonant relaxation between the disc stars. We 
shall also assume that there is no massive perturber in the 



disc. We may then focus on the effects of vector resonant re- 
laxation with the surrounding cluster. In the next sections 
we develop analytic machinery to provide a description of 
this interaction. 



3 LAPLACE-LAGRANGE THEORY FOR AN 
ISOLATED DISC 

Vector resonant relaxation affects the orientation of stellar 
orbits but not their semi-major axes or eccentricities. Since 
the relaxation time is much larger than the apsidal preces- 
sion time (|12p each orbit may be thought of as an axisym- 
metric planar annulus obtained from the Keplerian orbit by 
averaging over mean anomaly and argument of pericentre. 

We investigate the dynamics of the disc using a sim- 
ple model system. To construct this we shall assume that 
the stellar orbits in the disc are nearly coplanar (/ <^ 1) 
and nearly circular (e <C 1). In fact we shall make an even 
stronger assumption: that e,I <^ Aa/a where Aa is the 
typical difference in semi-major axis between a star and its 
nearest neighbor. This condition is satisfied in many plan- 
etary systems, including our own, but probably not in the 
Galactic-centre disc(s); nevertheless we shall argue below 
that it allows an analytic treatment of resonant relaxation 
that captures its most important features. 

With these approximations, we evaluate the Hamil- 
tonian of a system of A*' stars of masses, semi-major 
axes, inclinations, and nodes rrn, Ui, h, and Q.i, i = 
0, . . . , A'' — 1, to quadratic order in the eccentricity and incli- 
nati on. This is the classical L aplace-Lagrange secular the- 
ory iMurrav fc Dermottll 19991 ) . In this theory the inclination 
and node are decoupled from the eccentricity and apse. The 
evolution of the latter two elements is scalar resonant re- 
laxation, which we have argued is unimportant because of 
the relatively rapid apsidal precession induced by the stellar 
cusp. The evolution of the former is vector resonant relax- 
ation; the relevant Hamiltonian is 



JV-l iV-l 



miTTij (1) 
max(a„a,)"-''3/2(«-) 



El 

7j 



+ 



<h__qj_ 

H lit 



where onj = min(ai, a^)/ max(ai, a^), 

(1) f . _2 r cos Ode 



h\)>,{a) 



(32) 



(33) 



TT Jo (1 -2acos6'-f-a2)3/2 
is the Laplace coefficienlF°l. 7; = m\^'^ {GM,ai)^^'^ , and 

qi = "fili sin fli, pi = — 7i/i cos fli (36) 

^'^ In terms of the complete elliptic integrals of the first and sec- 
ond kind 



/■TT 

K{k) = / — 
Jo Wl 



Vl - sin^ ( 



m = / V 

Jo 



1 - A:2f 



(34) 



the Laplace coefficient satisfies 
.(1) 



V2(°) = -T^ ^ + - (1 ~ a^)K{a)] . (35) 

^' 7ra(l — a^)^ 
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are canonical coordinates and momenta. For some purposes 
it is useful to 'soften' the Laplace coefficient to suppress the 
singularity at q = 1, replacing the formula above by 



"3/2 



(a,6) = - 



cos 6 dO 



(37) 



TT J a (1 - 2a cos 6 + + e2)3/2 

where e is the dimensionless softening parameter and < 
e< 1. 

We use {x,y,z) to denote the standard Cartesian coor- 
dinates relative to which the node O and inclination / are 
measured. The total angular momentum is 

/ \ jv-i / sin/iSinSli \ 

I Ly I — mi\/GM,ai [ — sin/iCosSli 

\ cos li j 



i=0 
JV-l 



/ li sin Q,i 

= ^ m^VGM.a, [ -h cos ] + 0(/f 

1=0 



1 - 

^ 2 » 



(38) 



2 2-* 

It is straightforward to show that all three components of 
the total angular momentum are conserved (see discussion 
following eqs. I49H5T]) . 

The Hamiltonian can be rewritten as 



H{q, p) =p^Ap + q^Aq, 

where p'^ = (po, • • • , Piv-i), q^ = (go,--. 
N X N matrix A is defined by 



^ Gniirrijaij 

8 max(ai,aj)7i7j 
Grmmkaik ^(i) 



8 max(ai, ak)'^f 



{Olik) 



Note that 



JV-l 

^ A,;,7, = 0- 
3=0 



(39) 

, gjv_i) and the 
if i / j 

if i = j. (40) 
(41) 



Since A is symmetric, it can be diagonalized in the form 



OAO' 



(42) 



where O is orthogonal (O^ — O^^) and A is diagonal. The 
diagonal elements of A are the eigenvalues of A and the 
columns of O are the normalized eigenvectors (see Figures 
[SHl] below for properties of the eigenvalues and eigenvectors 
in simulated discs). 

Now consider a canonical transformation to new coordi- 
nates and momenta {Q,P) having the generating function 



S = P^O^q. 



(43) 



Then 



p=|^ = OP, Q=|^ = 0^q, q = OQ, (44) 



and 



H{P, Q) = P^AP + Q^AQ. 



(45) 



The inclinations of stars can be calculated from Q and P as 



2 I 2 

r2 1j + Pj 



7? 



7j 



7j 



(46) 



Since H is separable, +Qi is a constant of motion for all 
i — 0, . . . , N ~ 1. The equations of motion arc 



dH 



" dp 
which have the solution 



= 2A,P,, P, = 



dH 



= -2A,Q,, 



(47) 



(48) 



Q^{t) =Q^{Q) cos(2A,t) + P(0) sin(2A,t), 
Pi{t) = - Qi(0)sin(2Aii) + P(0)cos(2Ait). 

Thus Ai can be regarded as a frequency. 

Since 63^2 (ct) is positive, equation p2|) implies that 
H ^ 0. Thus H{q,p) is a positive semi-definite quadratic 
form defined by the matrix A, and thus all of the eigenval- 
ues Ai are non-negative. The lower bound H — is achieved 
if Pi/"/i = const, Qi/'yi = const independent of i, correspond- 
ing to a uniform tilt of the disc plane. Thus there is one zero 
eigenvalue, which we call Aq, and the rest are positive; in 
the discussions below we always order the eigenvectors by 
eigenvalue, so Aq = < Ai < - - - < Ajv-i- The eigenvector 
corresponding to zero eigenvalue is c(7o, . - . , 7Ar_i) where c 
is a constant (cf. eq. I41|) . Since the columns of O are nor- 
malized eigenvectors, Ojo = cjj with c = ±(5]] . 'y'j)^^^^ and 



Qo = ^ Ojoqj = c ^ JjQj = cL^, 

J=0 j=0 
N-1 N-1 

-Po = E '^^oP-' = ^ E ^jPj = ciy, (49) 

j=0 j=0 

where the last equalities follow from H38|) . Since Aq = 0, Po 
and Qo are constant, so the x and y components of the total 
angular momentum are conserved, as they must be- 
We define 



Z{q,p)^J2 



2 , 2\ 
Qi +Pi) 



(50) 



Following iLaskail l|2000l ). we call the 'angular-momen- 
tum deficit' since it represents the difference between the z 
component of the angular momentum of the actual system 
and the angular momentum that it would have if all of the 
stars were on coplanar circular orbits (eq. I38|) . Since the 
canonical transformation from {q,p) to {Q, P) is orthogonal 



ZiQ, p)^J2 (Q? + = + 



(51) 



Since + Pi is constant for all i, Z is also constant, which 
confirms that the z component of angular momentum is con- 
served, as it must be. 

Normal modes that have wavelengths that are short 
compared to the disc radius but long compared to the radial 
separation between stars (roughly, a 2> A 3> a/N where a 
is the disc radius and A'^ is the number of stars) approx- 
imatel y satisfy the WKB dispersion relation for b ending 



imatel y satisfy the WKd dispersion relation tor b endi 
waves l|Hunter fc Toomre|[l969l : iBinnev fc Tremaine| [200S 



(oj — mO,) 



2 , \2GE 
V + (27r) — 



(52) 



where v is the vertical frequency arising from an external 
potential and m is the azimuthal wavenumber. In our case 
m = 1, V = Q,, and <Si fi, so the dispersion relation 
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simplifies to 



OA ■ 



(53) 



The n' normal mode has approximately n nodes, and so the 
average wavelength of this mode is roughly 2(rniax — rmin)/fi- 
The longest waves have the lowest frequencies, and the fre- 
quency of mode n is proportional to uin oc n, consistent with 
the low- frequency behavior seen in Figure [3] below. 



4 INTERACTIONS BETWEEN THE DISC 

AND A SURROUNDING STELLAR SYSTEM 

Let us suppose that there is an external perturbation to the 
disc, which can be represented by generalized forces fqi, fpi 
that change the coordinate and momentum of ring i at a 
rate qi = fqi, pi = fpi. Since jiQi and 'jiPi are the x and y 
components of the angular momentum of star i (cf. eq. I38p 
we have 



(54) 



where T^iit), Tyiit) are the x and y components of the orbit- 
averaged external torque on star i. If these torques arise 
from a surrounding distribution of stars that is stationary 
and spherical on average, then {fqi{t)) — {fpi(t)) = for all 
i where (■) denotes an ensemble average over realizations of 
the surrounding stellar cluster. With the same assumptions 
we also have that {fqi{t)fpj{t')) = for all i and j. Since the 
number of stars in the surrounding cluster is large, fqi (t) and 
fpi{t) are Gaussian random processes. These can be uniquely 
characterized by the correlation coefficient 

T,,{t,t') = r,,i\t - = {fq^{t)fqAt')) = {fAt)UAt')). 

(55) 

where the first and last equality are valid if the forces are 
stationary and spherically symmetric on average. Note that 
Vij (t) need not be diagonal since nearby stars are expected 
to experience similar torques. The coherence time r of the 
external torques is defined so that 



Tij{M) ~ for |At|^r. 



(56) 



In the next two subsections, we discuss the general prop- 
erties of the growth of perturbations from an initially thin 
coplanar disc, then we construct a specific model for the 
torques and correlations. 



4.1 Response of the disc to external torques 

The equations of motion (|47l) in the canonical coordinates 
(Q, P) are modified to 

JV-l JV-l 

Qi = 2k,P, + Oj,fqj{t), P, = -2kiQ, + Oj,fpj{t). 

(57) 

Let us assume that the modes of the disc are non-degenerate, 
that is. A,; = Aj only if z = j (this is not an important re- 
striction for practical purposes). With the initial conditions 



QiiQ) ~ Pi{0) ~ 0, equations (|57|) have the solution 

Q^{t)^^Ojr / dti{fqj(U)cos[2K,{t-U)\ 

+ /pj(ti)sin[2A,(t-ti)]} 

= / dti{fp,(ti)cos[2A,{t~ti)] 

j=o Jo 

-/,j(ti)sin[2A,(t-ti)]}. (58) 

The mean squared value of the process over different real- 
izations of the perturbing torques is 



N-l 



= Y, OniOmi I dtl / dt2 
n,m=0 -^0 Jo 



m-Q " " 

X {{fqn{ti)fq,n{t2)) cos[2A,(t - ti)] cos[2A,(t - ta)] 
+ {/pn(ti)/pm(t2)) sin[2A,(f - ti)] sin[2A,(t - t2)]} 

OniOmi J dtl j dt2Tnm{t2-U)cOs[2Ki(t2-tl)\ 
n,m — Q 



= 2 ^ 0„«0™ / dt' {t-t')V^^{t')cos2Kt' (59) 

n,m=0 •'0 

where in the second-last line we have used the definition of 
Tnm(t) (ea. l55p and the trigonometric identity for the sum of 
cosines, and in the last line we have used Vnm{t) ~ r„,„(— t). 



{Q^{t)P^{t)) =0. (60) 

Similarly, we can calculate the cross-correlation coeffi- 
cient between two modes with Ai 7^ Aj, 



mt)QAt)) = {m)PAt)) 

N-l 



'mj / dtl / dt2 

n,m = ■'O JO 



n,m = 

X {(,f<7n(tl).f9m(fo)) C0s[2Ai(t - tl)] COs[2Aj(f - ti)] 

+ (/pn(ti)/pm(t2)) sin[2Ai(t - tl)] sin[2Aj(t - ta)]} 

= Yj OniOmj / dtl I dt2 
n,m=0 -^0 Jo 

Vr^m{t2 - tl) COs[2Ai{t ~ tl) - 2Aj(t - t2)] 

r-1 /* ,x ,t 



X 

AT 



— ^ ^ OniOmj ^ / dt rim{t ) 

n,m = Jo 

X [ sin(Aijt - 2Ait') + sin( Aijf + 2Ajt')] (61) 

where we have introduced A^j = Ai — Aj to simplify nota- 
tion. 

A similar calculation shows that 

{Q^{t)Pj{t)) = tan{A,jt){Q,{t)Qj{t)). (62) 

The mean squared inclinations at time t follow from 
equations ((46|, ([Mj, and (f6T|) . 

^-1 r) n 

u'w) = E liQAmit)) + {PAmm m 

k,i=o 
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tAi 



o 



harmonic 
oscillation 



independent 
random walk 



t > max(T, |Ay|-i) 



linear growth 



t < min(T, A, ^) 



correlated 
random walk 

r «: f < |A,j|-i 



temporally coherent -j^ temporally incoherent f^j- 
perturbations perturbations 

Figure 2. Evolutionary stages of the excitation of normal modes 
of an initially thin, flat stellar disc by stochastic external pertur- 
bations. On the horizontal axis, the age of the disc t is measured 
relative to t, the coherence time of the perturbations. The vertical 
axis shows t relative to the inverse frequencies A~ ^ of the normal 
modes of the disc. If t <S; r the perturbations are temporally co- 
herent, and the normal mode amplitudes initially grow linearly, 
then saturate and oscillate with angular frequency 2Ai . For t ^ t 
the perturbations are temporally incoherent, and the amplitudes 
of modes i and j undergo initially correlated, later independent, 
random walks (before and after time iA""*^! = |Ai — Ajl"-*^, respec- 
tively). The evolution in these regimes is examined separately in 
^iXn 11X21 and liXSl 



In the rest of this paper we explore the implications of 
equations (|57p - (|63p for the evolution of an initially planar 
disc excited by stochastic torques. First, we make general 
statements, then we examine applications to the Galactic- 
centre disc(s) in the next section. 

These equations contain three characteristic time- 
scales: (i) the age of the disc t, which determines the interval 
over which the external torques have acted on the initially 
thin and flat disc; (ii) the inverse frequencies A~^ of the 
normal modes of the disc; (iii) the coherence time r of the 
perturbations (eg. I56|) . 

Note that any perturbation that is initially zero, and 
then grows and decays smoothly with coherence time r (i.e., 
a perturbation whose temporal power spectrum contains 
only frequencies cannot excite short-period normal 

modes (A"'^ <C t) because the actions and hence amplitudes 
of these modes are adiabatic invariants (an analogy arises in 
the disruption of open clusters by molecular clouds: the tidal 
forces from the clouds are ineffective if the orbital period of 
the stars in the clust er is small compared to the passage 
time; see zeijll958f). 

The time evolution of the disc can be understood ana- 
lytically in the following limiting cases, shown schematically 
in Figure [2] 



4.1.1 Initial response, t <^ A^^ 

When the age of the disc is much less than the inverse fre- 
quency of a particular normal mode i, the factor cos 2 Ait' 
in the integrand of equation (|59|) is unity and we have 



(64) 



(Q?) = {Pf) =2 OmOrm / dt' (t - t')rnm(t') 

n,m = 



In particular, the zero-frequency (i = 0) normal mode de- 
scribing the mean orientation of the disc (cf. eq. I49p always 
grows according to equation (|64p . 

If the disc age is much smaller than A~^ for all modes 
i, then equation (|6ip yields 



N-i 

{Q^Q,} = {P'Pj} = 2 ^ OmOmj / dt' (t-t')rnm(t') 



and the orthogonal transformation to {qi,pi) yields 



(65) 



{q^ ) ^ (Pi) 



dt'{t-t')^^^{t'), ift<minA,"\ (66) 



(The same result would obtain if the disc mass were small 
enough that the collective effects from its self-gravity were 
negligible.) The inclination distribution is given by equation 
(HI, 

= A t dt'{t~t'){n,{0)nj{t')), t^inmA-\ (67) 
7, Jo 



4-1.2 Temporally coherent perturbations, t <^ t 

When the disc age t is much shorter than the coherence time 
r, the external forces are approximately constant over the 
lifetime of the disc, so T^jkit) ~ const, and equations (|59|l 
and (|61|) simplify to 



(0?) HP^} = c>„,o™r„ 

n,m — 
JV-1 



sin Ait 



„ , ■^-^ „ „ „ sin Ait sin A,t 

{Q^QJ) ={P^Pj) = Y OmO^jTum^^ -J- COS A,it. 



n,m — 



A, 



A. 



(68) 

Initially, the amplitude of each normal mode grows lin- 
early with time with a rate independent of Ai. Then 
after a saturation time fs,i = ^tt/A^ the amplitude 
reaches a maximum and begins to oscillate. At much larger 
times the time-averaged mean-square amplitude is (Qi) = 

2 y ^' nm ^niOmi^ nm / -^i • 

The time evolution for any particular realization of the 
perturbing forces can be obtained directly from equation 
(f58)l by specializing to constant fqj and fpj , 

JV-l 

= J2^ if"^ 2A^t - /„ (cos 2Ad - 1)] , 



3 = 
JV-1 



^'(*) = J2?r: [/" ™2A.t + ,f,, (cos2A,t - 1)] (69) 



2A. 



for i > (Ai / 0) ; otherwise for i = (Ao = 0) 

JV-1 JV-1 

Qo(t) Ojo/,., Po{t) = tY,0,ofpj. (70) 

i=o j=o 
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4.. 1.3 Temporally incoherent perturbations, t ^ r. 



or wires of uniform density: 



When the disc age is much larger than the coherence time, 
then since Tnm{t') ~ for t'^r, the upper hmit of the in- 
tegration in equations (|59l) and (|61|) can be truncated at r, 
so 

(0? W) HPht)) = ga - h,, if t;>r, (71) 

COS " t 

{QiifilQiii)) = . [^ij cos Aijt + Bij smArjt] , 
ij cos Aijt + Bij sin Aijt] , 



where the constants Aij, Bij, gi, and hi are 

JV-l 

A,, = ^ 0„«0,„j / dt' r„^(t')[sm{2Ajt')-sm{2A,t')] 

n,m=0 ''0 
JV-l 

-Bij = OniOrnj / dt' (t') [cOs(2Ai ) + COS (2Aj t')] 



/o 

JV-l 



Qi ^Bii, hi =2 ^ OniOmi / dt' t'rnmit')cOs2Ait' . 
n,m-0 

(72) 

For t > max(r,ft,/30, = {Pi{t)) ^ 9^t■, 

moreover {Qi{t)Pi{t)) = (eq. [60]), and {Qi{t)Qj{t)) and 



{Q^{t)Pj{t)) are bounded by iy'A?^. + B^ / jA^^ j, so after 

t ^ max \t, |Aij|^^(^ij + Bfj^^^/gi^ , the cross-correlation 

becomes negligible compared to {Ql{t)). In other words, at 
large times Qiit) and Pi{t) undergo random walks. The ran- 
dom walks of various modes i and j are first correlated then 
gradually become independent after the disc lifetime be- 
comes much larger than the inverse relative normal mode 
frequencies |Aij|~"'. 

In our case none of these simple limits applies: the co- 
herence time of the perturbations from the cluster is the vec- 
tor resonant relaxation time-scale, and this can be shorter 
or longer than the disc age depending on the radius and the 
effective mass m2 in the cluster (cf Fig. Moreover the 
set of inverse frequencies includes values that are both 
shorter and longer than the disc age (Fig. [3]). We compare 
these time-scales numerically for Monte Carlo simulations 
of the stellar disc in ij5] below. 

Finally, the configuration of an isolated disc in thermal 
equilibrium under vector resonant relaxation is discussed in 
the Appendix. 



4.2 Distribution of torques 

We now construct a model for the perturbing torques 
{Tjji, Tyi}. These are the sum of the torques from all of the 
stars in the cluster on the disc star labelled by i, in a circu- 
lar orbit near the z = Q plane; the torques are averaged over 
the orbit of both the cluster star and the disc star. To sim- 
plify the calculation at modest cost in realism, we assume 
that the cluster stars are also on circular orbits. Then the 
total torque on disc stars on circular orbits can be calculated 
from the gravitational torque acting between circular rings 



Ti = 'YY^KipiP2t{ni3-ni)niXni3, (73) 

p 1=1 

where i and /3 label disc and cluster stars, respectively, Li 
or Lp is the angular momentum vector of star i or /3, and 
np = Lp/\Lp\ is the unit vector aligned with the angular 
momentum, which is related to the node Q.^ and inclination 
//3 by n/3 = (sin sin — sin /;3 cos r2/3, cos J/j). We also 
have 



Ki^i = Gmim^[P2pXQ)] 



2 'i/9< 
^2^+1 ■ 



P2l(Q) 



(-i)'^r(2^ + i) 



22£r2(£+l) ■ 

(74) 

Here and rp axe the mass and orbital radius of star /? (re- 
call that both the disc and cluster stars are assumed to be on 
circular orbits), ri/3< = min(ri,r/3), ripy = max(ri,r;3), and 
P2e(x) and P2i(x) are the Legendre polynomial and its first 
derivativ^HI- The sum over I in equation (|73|) is absolutely 
convergent unless = rp, The analytical formula (|29|) for 
the precession rate when r^/ry ^ 1 is the corresponding 
limiting case of equation (|73[) . 

In the discussions below, we neglect the back-reaction 
of the disc stars on the orbits of stars in the old cluster. This 
approximation should be valid so long as the total angular 
momentum contained in the fiuctuating non-spherical com- 
ponent of the star cluster exceeds the angular momentum 
in the disc. The former is roug hly N^''^m2{GM.rY''^ where 
7712 = {m?) / {m) is the effective mass and A'" = M{r)/m2 
is the effective number of stars in the cluster at radius r; 
the latter is M^isc{GM,r)^^^ . Thus the neglect of the back- 
reaction should be valid so long as [Af(r)m2]^^'^>Mdisc. This 
is the same as the criterion that nodal precession is domi- 
nated by the stochastic field from the old cluster rather than 
the disc, as discussed following equation (|26[) . 

Assuming that the disc is fiat (i.e., zero inclination for 
all disc stars), equation (|73|l simplifies to 



Til' 



^ 5Z ^^f"P^t (cos /a ) 



cos Q.fj 
sin Q,fi 



(75) 



4-2.1 The torques as a Gaussian random process 

In a spherical cluster, or the equatorial plane of an axisym- 
metric cluster, the torques described by equation (|75p are 
the sum of a large number of independent random variables 
with zero mean. Hence, according to the central limit theo- 
rem, the probability distribution of the torques T^i and Tyi 
is Gaussian, with zero mean and covariance or correlation 
function (T,»(0)T,j(t)) = {Ty^{0)Ty,{t)) = 7«7jT»j(t). From 
equation (I75p it is also clear that (TxiTyj^ = in a spherical 
cluster, for all i and j. 

A Gaussian random process is fully characterised by 
its mean and correlation function. Thus, when generating a 
Monte Carlo distribution of cluster stars, the distribution of 
torques is the same for all mass functions having the same 
effective mass m2 so long as the total mass distribution M{r) 



In terms of the associated Legendre function P^{x), 
P2g(cosI) = — P2^(cos/)/ sin/. 
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is fixed. Tfius we may assume witliout loss of generality that 
all the cluster stars have mass m/3 = m2. 

Next we examine the spatial correlation function of clus- 
ter torques (r„(0)r,j (0)) 

— lilj^-iji^)- We are mostly in- 
terested here in the correlation function within a thin, flat 
disc. Then we may use equation (|75p for the torques, which 
yields 

{T,,(0)T,,(0)) 



= ( ^ Kif3lKjfii,P2t{cOsIfj)P2t>{cOsIl3)siV? 0,13 



min(ri, rj)/ max(ri 



„2* 



„2^ 



'i/3< I jl3< 



/3 



where ke is a constant defined as 



2C+1 2e+i ' 



(76) 



_ 2£il + 2£) r^(£ + l/2) 

7r(l + 4^) r2(^-|-l) • ^ 

Here we have used the orthogonality of the Legendre func- 
tions, (P27(cos/^)P2" (cos/;,)) = for ^ / £'. 

We have separated ke from the factor [P2«(0)]^ in equa- 
tion (|76p because ke is practically independent of £; as £ 
varies from 1 to oo, ke varies only from 3/10 = 0.3 to 1/n — 
0.3183. Therefore we can estimate the correlation function 
to reasonable accuracy by replacing ke with a constant k, 
0.3 ^ fc < 0.32. This leads to a power series in £ which is 
reminiscent of the power series defining the complete ellip- 
tic integral of the first kind, K{x) = ^tt J]^^JP2n(0)]2a;^". 
Therefore 



(r„(o)T,,(o)) 

= G^kniimj 



(78) 



where atfj = rip^/rip^ = min(ri, r/j)/ max(ri, r^). We may 
express this result in terms of the density of the cluster p(r) 



(eq. [7| and the effective mass m2 

oc 

dr r^p{r) 



max(r, n) max(r, rj) 



'-K{aiaj 



(79) 



where a; = min(r, r^)/ max(r, r;). 

Equations (|78|l - (|79|) define the probability density of 
torques acting on the disc stars at each instant. They 
show that the correlation function for disc stars at radii 
Vi and rj is determined mainly by cluster stars in the re- 
gion min(ri, rj)<r< max(ri, rj). In particular, at large radii 
p(r) cx r"^'^^ (eq. [7| so for r 3> max(ri,rj) the integrand 
declines as rfr|r~^'^^, so the contribution to the correla- 
tion function from this region is negligible. At small radii. 



p{r) 



so for r ^ min(ri,rj) the integrand scales as 



^ ' /{fi'i'j), impl ying that the co ntribution from small radii 
is also negligible. iGiirkan fc Hop man ( 2007) reach a similar 
conclusion. 

A simple fitting formula for equation (|79[) . using the 
mass distribution in the Galactic centre (eq. |8]), is 



{T,«(0)T,,(0)) ~CT 



m2p{rij)r% „ 
■'-UL,Q.{r,)n(rj)a, 



(80) 



where ai 



miriD,{ri) is the angular momentum, and ct and k are fit 
parameters. The best fit values, which approximate equation 
H79|) to better than 20% over radii Oij < 4 and fij < 4pc, are 
Ct ~ 0.77 and k — 1.8. Thus, the correlation of torques on 
different stars, (T^iT^j) /[{T^iy^^{T^jy^'^], is less than 14% 
for Ui ^ 3aj. Furthermore, the distribution of torques is a 
smooth function of radii with a peak at ai — aj. 



5 MONTE CARLO SIMULATIONS 

We construct a fiat, razor-thin disc of TV stars on circular or- 
bits. The disc is assumed to lie initially in the reference plane 
so li = Qi/yi = Pi/^i = 0, for all i = 0, . . . , A'' — 1. The semi- 
major axes are chosen randomly from the surface densit y 
distribution ([1} with exponent 5 = 1.4 (jBartko et al.|[201ol ). 
between inner and outer radii of 0.04 and 0.6 pc (roughly 
larcsec to 15.5arcsec). The stellar masses are chosen from 
the mass function (|14p with a — —0.45, mmax = 30 Mq and 
nimin = 1 M© ; the minimum mass is arbitrary but this choice 
has almost no influence on our results, and the maximum is 
the most massive star that can survi ve for the 6 Myr age of 
the disc(s) l|Leieune fc Sch aerer'200lh. We set the total num- 
ber of stars to be 500 in this mass range, implying that ~ 120 
stars have masses M 20 Mq , consistent with observations 
- 90 massive WR/O stars h ave been observed with ~ 75% 
spectroscopic completeness (|Bartko et al.l 120091') . and thes e 
typically have masses m > 20 Mq ( Paumard et al.l |2006| ) . 



The corresponding disc mass is 6.3 X lO'^ Mg 



5.1 Normal modes 

As discussed above, the evolution of the disc is most easily 
described in terms of its normal modes. We generate 1000 
Monte Carlo realizations of the disc (stellar masses and semi- 
major axes) . For each realization, we calculate the matrices 
A and O fegs. 1401 and I42p . as well as the eigenvectors (the 
columns of O). As explained in ii4.1l the evolution of the disc 
is determined by the relation between three characteristic 
time-scales: the inverse normal mode frequencies A~^, the 
coherence time r of torques from the spherical cluster, and 
the age of the disc t. 

Figure shows the distribution of the saturation time 
ts,i ~ ^7r/Ai, the characteristic time-scale at which the am- 
plitude of a mode subjected to a fixed torque stops growing 
and begins to oscillate (see ^4.1.2^ . The bottom (blue) curves 
show the mean distribution and the 95% confidence interval 
for our standard disc model. The figure shows that all but a 
few modes satisfy ts.i < 6 Myr, and so are already saturated 
at the current age of the disc. 

Some sample normal modes are shown in Figure [4] The 
height of the orbit of star j above the reference plane at az- 
imuth (f> is Zj — rjlj sm((f)—Qj) = —'y^^rj{qj cos 4>+Pj sinrf))- 
Thus if only normal mode i is present, with amplitude 
{Qi,Pi), the height is Zj = —rj{Qj coscp + Pj sm(l))Oji/'yj. 
The figure plots Oji/jj which is proportional to the frac- 
tional height Zj/rj at fixed azimuth. The left panel shows a 
single realization of our standard disc model and the centre 
panel shows an average over 1000 realizations. The modes 
are ordered by increasing frequency; mode i — (the black 
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Figure 4. Normal modes of the stellar disc i = 0, 1, 2, 3, 4, and 16 (ordered by increasing frequency) for a randomly chosen realization 
(left), and averaged over 1000 realizations of the disc (middle). The vertical axis is proportional to the fractional height of the mode 
or the inclination angle at fixed azimuth. These modes of the stellar disc oscillate independently in the absence of the cluster. Right: 
Probability distribution of the saturation times for the same modes. A vertical green line shows the disc age t = 6 Myr. The models 
contain 500 disc stars with masses between 1 and 30 and radii between 1 arcsec and 30 arcsec. 
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Figure 5. Driving modes of the spherical cluster, i.e., the eigenvectors of the spatial correlation function of the torques on the disc, 
(TxiTxj). The figure shows modes i = 0,1,2,3,4, and 16 (ordered by decreasing significance) for a randomly chosen realization (left), 
and averaged over 1000 realizations of the disc (middle). The models contain the same stars as in Fig. |4l Right: Probability distribution 
of the mean squared amplitude for the same modes. 




7r/(2Aa [yr] 

Figure 3. Histogram of saturation time-scales Ivr/Ai for a disc of 
stars distributed between 1 arcsec and 15.5 arcsec (dashed blue) or 
extrapolated to 525 arcsec (solid red) with surface density S(r) oc 
J,— 1.4 rpj^g stellar masses are chosen from equation (I14I I with a = 
—0.45, mmax = 30 Mq, and mmin = IMq; the total number of 
stars is 500 and 5000 in the two cases. Dotted lines show the 95% 
confidence interval. The peak of the distribution corresponds to 
eigenmodes with average wavelengths comparable to the average 
distance between neighbouring stars. The age of the disc(s) in the 
Galactic center is marked by a vertical green line. 



horizontal line in the figure) has frequency Ao = (eg. I49p . 
and corresponds to a uniform tilt of the disc. 

The low-frequency modes are rather smooth, and their 
shapes are approximately the same in different realizations 
of the disc - the mode shapes are not sensitive to the 
large random variations in the masses of individual stars 
(1 Mq ^ rrij ^ 30 Mq). The long- wavelength modes are well 
represented by the average waveform, while shorter wave- 
length modes in a single realization deviate more strongly 
from the mean (compare the orange lines, i = 16, in the 
left and centre panels of Figure S} . The right panel shows 
the probability distribution for the saturation times ^iv/Ai 
for various realizations of the disc. The probability distri- 
bution is sharply peaked for Ai (for each fixed i) with a 
FWHM of about 20%. Therefore, Ai can be predicted from 
the surface-density distribution without knowing the loca- 
tions and masses of disc stars. The normal modes i = 0, 1, 2 
are still growing {ts,i > 6 Myr), i = 3 is just around satura- 
tion, and all other modes are saturated. 

These considerations and the theory presented in H4.1.2I 
allow us to make general remarks on the expected warping 
of the disc. Recall that modes in the oscillating phase - 
with either the coherence time or the disc age less than the 
saturation time ts,i — ^tt/ Ai - are suppressed in amplitude 
relative to lower frequency modes by A~^ (see eq. [68] and 
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i]5.3l below'). Therefore if the external torques on the disc 
stars do not depend too strongly on radius, the shape of the 
disc is expected to be dominated by long-wavelength modes. 

There are presently few observational constraints on the 
maximum radial extent of the disc, and its evolution would 
be different if it extended to larger radii. As an example, the 
top (red) curves in Figure[3]show the histogram of saturation 
times for a hypothetical disc extrapolated to 20.4 pc with the 
same power-law surface density distribution. In this case, 
there are many more modes still in the growing phase at 
the current age of 6 Myr. However, the characteristic wave- 
lengths of these additional modes are larger than 0.6 pc, and 
the saturation timescales of the smaller wavelength modes 

are not very different from the case in which the disc is trun- — 
cated at 0.6 pc. Therefore we expect that our predictions of 
the shape of the warped disc at radii J50.5 pc are robust, and 
independent of whether or not there is an outer disc. 



5.2 Spherical cluster 

We model the torques on the disc from the star cluster as 
follows: (i) As described above, since the torque distribu- 
tion is a Gaussian random process, we may assume that all 
cluster stars have the effective mass mi, as defined follow- 
ing equation (|13p : we typically assume mi = 10 M©. (ii) As 
described above, for the sake of simplicity we assume that 
the cluster stars are on circular orbits, and generate a distri- 
bution of orbital radii consistent with the mass distribution 
M(r) in equation ((Sjl out to a radius of 2pc, sufficiently far 
outside the disc(s) that the perturbations from larger radii 
should be small, (iii) We assign nodes f2a uniformly ran- 
dom between and 27r and inclinations //3 so that cos 7/3 is 
uniformly random between —1 and consistent with an 
isotropic distribution of orbits, (iv) We expand the torque 
(|73p to order 21 = 10. (iv) We assume a disc age of 6 Myr. 

Alternatively, random realizations of cluster torques 
may be generated by sampling an A^-dimensional Gaussian 
distribution (where N is the number of disc stars) with a co- 
variance matrix {T^iT^j) given by equations (|79p - (l80|l . We 
can decompose the torques into independent driving modes, 
corresponding to the eigenvectors of (TxiTxj). These driving 
modes, shown in Figure O are a property of the spherical 
cluster of old stars and hence are distinct from the nor- 
mal modes of the disc shown in Figure |4l We find that the 
corresponding eigenvalues, which correspond to the mean 
squared amplitude of the particular driving mode, span a 
vast range, some 6 orders of magnitude. Remarkably, the 
five longest wavelength driving modes typically contribute 
90% of the total torque. Therefore, the cluster excites the 
disc predominantly through a few long-wavelength driving 
modes. 

As described at the end of i]4.1.3l the coherence time 
of the torques from the cluster can be shorter than the disc 
age depending on the radius. The two are equal at a radius 
around Tt ~ 0.3 pc(m2/10 Mq)"'^^. Nevertheless, to avoid 
excessive complication - both conceptual and numerical - in 
our simulations we shall assume that the coherence time is 
long compared to the disc age, so that the perturbing forces 
fqi{t), fpi{t) can be regarded as time-independent and the 
formulae in tj4.1.2l can be used to calculate the disc evolution. 
This simplification has a smaller impact on the evolution of 
long wavelength modes A>rr, while it may be important 
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Figure 6. Two orthogonal cross-sections of a disc subjected to 
vector resonant relaxation. The parameters of the disc and the 
stellar cluster in which it is embedded are described in ^ The 
disc was initially horizontal. The calculation is based on the as- 
sumption that the inclinations are small and the results presented 
here are large enough that this approximation is not quantita- 
tively accurate. Only stars with mass > 20 Mq are shown but the 
stars of lower mass have the same thin, warped distribution. The 
simulation assumes that the disc age t = G Myr and the effective 
mass in the stellar cluster is mi = 10 Mq; the inclinations scale 

1/2 
as mj . 

for short wavelength modes close to the center. To account 
properly for the temporal evolution of the torques from the 
cluster would require solving simultaneously for the secular 
evolution of both the cluster and the disc stars (see ^ . 

5.3 Warped disc 

The results of a typical simulation are shown in Figure 
[S] which plots two orthogonal cross sections through the 
disc. The disc remains thin, but exhibits a substantial warp 
with a structure dominated by the zero-frequency tilt mode 
and the three longest wavelength normal modes of the disc 
j — 1, 2, 3. We conclude that vector resonant relaxation can 
excite a coherent warp in an initially thin, flat disc. 

In two-body relaxation the perturbing forces from a sur- 
rounding cluster of stars usually thicken a disc, rather than 
warping it. For example, the observed thickness of the Galac- 
tic disc in the solar neighborhood has been used to constrain 
the effective mass of t he objects comprising th e Galaxy's 
dark-matter halo (e.g.. iLacev fc OstrikeJ [l985l ) . Warping, 
rather than thickening, occurs in this case for two main rea- 
sons: 

(i) Vector resonant relaxation arises through the torque 
from a stellar orbit after averaging over mean anomaly and 
argument of pericentre, and this averaged torque has much 
less small-scale power than the forces exerted in the passage 
of a nearby star (cf. Fig. [5]). More quantitatively, consider 
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Figure 7. Identical to Figure (6] except that the gravitational 
interactions among the disc stars have been turned off. The disc- 
star inclinations are much more irregular at small radii. 



two nearby stars on circular orbits with semi-major axes a 
and a -I- Aa, and relative inclination I<Aa/a; torques from 
stars at larger relative inclinations are smaller and can be 
neglected for this argument. The time-averaged or secular 
torque between the two orbits is T ~ Gm'^/Aa. For a radial 
density profile n ~ nQr~'' , the number of stars in this semi- 
major axis and inclination range is AA^ ~ noa^~^ Aa ~ 
nga '{Aa)^. The total stochastic torque generated by these 
stars is then {AN)'^^^T oc nl^^Gm^{Aa)'^^^ . This is an in- 
creasing function of Aa, showing that the secular torque on 
the disc is dominated by large-scale fluctuations. A nonzero 
eccentricity for either the cluster or disk stars further re- 
duces the small-scale power. 

(ii) As suggested by equations ((58]), the excitation of 
normal modes in the disc is reduced if the coherence time 
of the torque exceeds the inverse frequency of the normal 
mode, and the coherence time for vector resonant relaxation 
is longer than the inverse frequencies associated with small- 
scale normal modes in the disc. To illustrate the importance 
of this effect, Figure[7]shows the results of a simulation iden- 
tical to that leading to Figure [51 except that the masses of 
the disc stars have been reduced to nearly zero. Reducing the 
disc-star masses decreases all oscillation frequencies Ai, since 
these are proportional to mass in Laplace-Lagrange theory. 
In the limit of near-zero mass, the interaction between the 
disc stars becomes negligible as described in Sec. 14.1.11 and 
{QiQj) becomes independent of A;. Figure [7] shows that in 
this limit the warp becomes much more irregular and less 
spatially coherent, because both large-scale and small-scale 
modes are excited. 

Figure |8] shows the normal mode power and cross- 
correlation as a function of the average wavelength, Ai = 
2()"max - ?-min)/i. The j/-axis shows (Qi+kQi-k) for k = 0, 
1, and 2, where Qi is the amplitude of the i"" normal mode 
(see Fig. [4]). The circles show the mean of and the shaded 
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Figure 8. The typical amplitudes of normal modes (arbitrary 
scale), as a function of their average wavelength. The circles mark 
the mean of the squared amplitudes of normal modes for a random 
realization of the disc for each i, and the 95% confidence interval 
is shaded. Plusses and crosses show the cross-correlation between 
various normal modes for comparison. In this case, the distribu- 
tion is symmetric around 0, and the symbols are placed at the 
boundary of the 95% confidence interval. The modes i = 1,2, 3, 4, 
and 16 shown in Fig.|4]are labelled. The torques are assumed to be 
constant during the 6 Myr lifetime of the disc with an amplitude 
corresponding to the RMS torque (eg. I80|l. The high-frequency, 
short-wavelength normal modes are suppressed. 



region shows the 95% confidence interval for a random re- 
alization of the disc. In contrast, the cross correlation of 
normal modes has zero mean, for these elements the sym- 
bols correspond to the upper boundary of the 95% interval. 
The calculations are based on equations (|68p and (|80|) . The 
long-wavelength modes dominate the disc, and modes with 
wavelength ^larcsec are suppressed by some five orders of 
magnitude. Note that the predicted amplitude of the longest 
wavelength modes {i 5C 6) exhibits only a small scatter be- 
tween various realizations of the disc. Although the different 
normal modes are not independent, the cross correlation be- 
tween modes with increasingly different average wavelengths 
are clearly more and more suppressed. The RMS warping of 
the disc then follows from eq. (|63|) . Figure [8] shows that a 
coherent torque tilts and warps the disc in a way that its 
shape is reminiscent of only a few long wavelength modes. 



6 DISCUSSION 

The massive young stars observed in the Galactic centre at 
radii ~ 0.05-0.5 pc from the central black hole may have 
formed in a thin, flat disc. Using the best available esti- 
mates for the mass, age, and other properties of the disc 
and the cluster of old stars in which it is embedded, we 
have shown that the disc naturally and inevitably develops 
a strong warp through vector resonant relaxation with the 
cluster. We suggest that this mechanism explains the large 
(~ 60°) warp that is observed in the disc. 

We have modelled the gravitational torques among the 
disc stars and between the disc stars and the surround- 
ing old cluster. Our models are based on Laplace-Lagrange 
theory, in which the secular gravitational interactions of 
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the disc stars are mathematically equivalent to the interac- 
tions in a system of masses connected by springs. Laplace- 
Lagrange theory assumes that the relative inclinations are 
small, Ajfc = \a-jlj — afe/fc|/|aj — a^l <C 1. For the simu- 
lation shown in Figure [6l the maximum value of Aj^ was 
2.8, indicating that Laplace-Lagrange theory has been ex- 
tended beyond its domain of applicability but should still 
be qualitatively correct. 

In Laplace-Lagrange theory, the disc dynamics can be 
decomposed into normal modes. The eigenfrequencies of the 
normal modes in our models span a vast range of time- 
scales (cf. Fig. (Sj. The high-frequency normal modes are 
an unphysical consequence of the assumption that the disc 
orbits are precisely circular and that the inclinations are 
small compared to the separation between stars; in fact the 
highest frequencies correspond to modes involving only two 
adjacent stars. However, our results should not be affected 
by this shortcoming since the high-frequency modes are not 
excited efficiently by torques from the cluster stars (cf. Fig. 
[HI). 

The frequencies depend monotonically on the wave- 
length of the mode: oscillation times tt/A of about 13Myr 
correspond to wavelengths of ~ 0.4 pc (lOarcsec) (cf. blue 
curve in Fig. 4); faster or slower modes have smaller and 
longer wavelengths, respectively. The frequencies and the 
shapes of the slow normal modes are insensitive to the de- 
tails of the stellar mass function so long as the large-scale 
surface density of the disc is fixed. 

The growth of a given mode saturates after about a 
half oscillation period. The low-frequency, long-wavelength 
modes are excited to much larger amplitudes than high- 
frequency, short-wavelength modes, both because the secu- 
lar torques couple more strongly to long-wavelength modes, 
and because the low-frequency modes saturate after a longer 
time (an equivalent statement is that the amplitudes of the 
high-frequency modes are adiabatic invariants). 

The most important free parameter in our models is 
the effective mass ni2 = {m^)/{m) of the stars and other 
irregularities (e.g., clusters, gas clouds, etc.) in the cluster; 
the simulation in Figure [6] assumes m2 — 10 Mq , which is 
likely to be an underestimate. A larger value of m2 increases 

1/2 

the warp amplitude proportional to . 

For simplicity our simulations are done using time- 
independent torques from a cluster. For a spherical cluster, 
this approximation is valid if the disc age is less than the 
vector resonant relaxation time, since this is approximately 
the same as the coherence time. A rough measure of the 
validity of this approximation is given by Figure [71 which 
shows the evolution of the orientations of the stars in a disc 
of zero mass: since the disc self-gravity has been turned off in 
this figure, the stars behave in the same manner as random 
stars in the spherical cluster, so the approximation that the 
perturbing forces from the cluster are constant is valid if the 
inclinations of the stars in this Figure are small. Evidently 
the approximation of constant forces is good for the outer 
parts of the cluster, where the relaxation time is relatively 
long, and suspect in the inner parts. 

Both of our major approximations - Laplace-Lagrange 
theory and constant perturbing forces - are more accurate 
if the effective mass is smaller than the assumed value m2 = 
10 M0 and worse if it is larger. 

We also assumed that the cluster is spherical on av- 



erage (i.e., apart from Poisson fluctuations due to individ- 
ual stars), that the orientations of cluster star orbits are 
uncorrelated, that there is no back-reaction from the disc 
on the cluster, that the cluster stars are on circular orbits, 
and that the cluster mass profile M{r) (eq. [8]) can be de- 
termined from the luminosity distribution L{r) assuming 
constant mass-to-light ratio. Under these approximations, 
the fractional deviation from sphericity at radius r>0.2 pc 
is roughly y/m2/M{r) = 0.009 V"i2/10 M0(r/O.l pc)-° '^^ 
The deviation from sphericity may be larger if the clus ter is 
flattened due to rotation; however, iTrippe et al.l ()2008 l) esti- 
mate the rotation speed of the cluster to be v^otir) = (3.6 ± 
0.8) km s~^(r/0.1pc) within Ipc. The rotational flattening 
of the cluster is then roughly Wrot/""^ ~ 2 x 10^*(r/0.1 pc)^, 
smaller than the stochastic flattening for r<0.3pc. The flat- 
tening of the overall mass distribution due to the disc is 
comparable to the stochastic flattening (see discussion fol- 
lowing eq. I26|) so the neglect of this contribution probably 
does not seriously invalidate our results. The justification for 
neglecting the back-reaction of the disc is described at the 
end of i]4.2l The assumption that the orientations of cluster 
star orbits are uncorrelated is harder to justify, and to do so 
will require a self-consistent simulation of both the disc and 
cluster stars (see the discussion at the end of this section). 
The assumption that the stars in the cluster are on circular 
orbits is unrealistic - the scalar resonant relaxation time is 
much less than the age of the Galaxy so we expect the ec- 
centricity distribution to be dn oc e de. This defect will be 
rem oved in future calculatio ns. 

iBregman fc Alexander! (|2009l ') have suggested that 
warps in accretion discs surrounding black holes in the cen- 
tres of galaxies - in particular, the ~ 10° warp in the maser 
disc in NGC 4258 - may arise from resonant relaxation 
with the surrounding stellar cluster. In accretion discs the 
warp dynamics may also be affected by gas pressure, viscous 
dissipatio n, and radial transport of mass and angular mo- 
mentum (|Pringlelll993 : iLodato fc Pringlell2007l ). as well as 
opacity and rad i ation pressure in radiatively efficient discs 
(jPettersonlligrrl : iMalonev. Begelman. fc Pringlelll996l '). 

The treatment in the present paper does not address 
several important issues: 

• What are the characteristics of the disc after vector 
resonant relaxation excites a large warp - in particular, large 
enough that the Laplace-Lagrange treatment is invalid? 

• Can vector resonant relaxation create structures that 
resemble the second, 'counter-clockwise' disc seen at the 
Galactic centre? A mechanism to form both the clockwise 
and counter-clockwise discs from a single thin, fiat precursor 
would explain why stars in both discs have the same age. 

• The WR/O stars in the disc(s) are much more massive 
than the old cluster stars. What is the analog of dynamical 
friction for vector resonant relaxation and what role does it 
play in determining the disc structure? 

These questions can be addressed most effectively by 
'N-ring' simulations that follow the evolution of a set of 
A'' axisymmetric rings, each representing the smeared-out 
mass density in a Keplerian orbit after averaging over mean 
anomaly and argument of pericentre. Each ring exerts a 
torque on every other ring, and the simulation follows the 
precession of the rings in response to these torques. We have 
written a code to carry out these simulations and are cur- 
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rently using it to follow the evolution of stellar systems that 
resemble the central parsec of the Galaxy. 
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APPENDIX A: THE MICROCANONICAL ENSEMBLE IN LAPLACE-LAGRANGE THEORY 

We have argued in !j2]that vector resonant relaxation among the stars of the Galactic-centre disc(s) is unimportant. However, 
there is a wide range of disc ages and environments for which internal vector resonant relaxation may be the dominant 
process of dynamical relaxation. A particularly simple case is an isolated disc with age much larger than the vector resonant 
relaxation time. Isolated discs conserve total energy and angular momentum; moreover, if vector resonant relaxation is the 
only important dynamical process the semi-major axis and eccentricity of each star is conserved. According to the usual 
principles of statistical mechanics, the probability distribution of such discs (the microcanonical ensemble) should be uniform 
in the manifold of phase space that is determined by these conserved quantities. 

For sufficiently small inclinations, the Hamiltonian of the disc is given by the Laplace-Lagrange form (|32p . Since this 
Hamiltonian is quadratic in the coordinates and momenta, the equations of motion are linear and therefore integrable. Thus in 
principle an isolated disc with sufficiently small eccentricities and inclinations exhibits no resonant relaxation. Our assumption 
is that higher-order terms that we have neglected in the Hamiltonian lead to relaxation to an equilibrium state that can be 
approximately described using the statistical mechanics of the quadratic Hamiltonian - just as occasional collisions in an ideal 
gas lead to a thermal equilibrium that can be described using the quadratic Hamiltonian H — ^p^ /m. 

The Hamiltonian H and the angular-momentum deficit Z (eq. I50p are both conserved. An important combination of 
these is 

r(q,P) = ^77 r = ^ ; — ^ • (Al) 

Z{q,P) p^p + q^q 



Since F has units of inverse time we call it the frequency parameter of the disown. 
We can re-write the last of these expressions as 

which shows that F can only span a limited range of values, from Amin to Amax where these are the minimum and maximum 
eigenvalues of A. 

We may also assume that the conserved x and y components of the angular momentum L^, Ly are both zero, since this 
can be ensured by choosing the z-axis of the coordinate system to be parallel to the total angular-momentum vector. Then 
from equation (|49|) Po ~ Qo ~ 0. 

To analyse the properties of the microcanonical ensemble, we first find the density of states ujm{E,C), defined so that 
ljn{E, C)dEdC is the volume in 2{N — l)-dimensional phase space in which the energy and angular-momentum deficit are in 
the ranges {E, E + dE) and (C, C + dC). ThuiQ 



ujn 



We use the identity 



Then 



„ N-l 

{E,C)= / Y[dP^dQr5[E ^ H{P,Q)]S[C ~ Z{P,Q)] 

„ N-l 

= Y[dP^dQ,5{E-P'^^P + Q^^Q)5{C~P'^P-Q^Q). (A3) 
5{x) = — / die'". (A4) 



ljn{E,C)=-— YldP^dQ^ dt dt' e-yir,[tE + t'C -t{P'^\P + Q^NQ)-t'{P'^P + Q^Q)\ 

J i — 1 —zoo J —too 

1 r pioo rioo ^ N — 1 

= UdPdQj die"" dt'e"^exp - ^(P." + Q?)(tA.+t') • (A5) 

■' i = l J -ice J-ica 

The next step is to exchange the order of integration, but this cannot be done immediately as J dPdQ exp 

is not absolutely convergent for imaginary t and t' . However, using the identity 1 = exp(MC) exp[— u5]]^~^(Pi'^ -f QI)] (recall 
Po = Qo = 0) we can rewrite the double integral as 

un{E,C)^-^ [ Y[ dPrdQ, Hdie"' T'df' e''^ exp [ - Y.{Pf + Ql)(tK, + t' + u)V (A6) 



Note that F used here is not related to Tnm{t) denoting the correlation function in the main text. 

The system we are studying resembles the celebrated spherical model for a ferromagnet llBerUn fc Kadll952l l and much of our analysis 
is borrowed from the literature on that problem. 
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and the integral over dPidQi is now absolutely convergent for u > 0. The order of integration can now be exchanged so that 
the 2(7V — l)-dimensional integral over phase space can be done: 

N — Z ricG r-u-\-ioo N—1 N — 3 fioa ru-^-ioc ^ N — 1 

ujn{E,C) = — / dte^ dse"'^ Ylis + tAi)-^ = — / dt ds exp + sC - ^ log(s + tAi) , 

J —ioo J u — ioo i — l —ioo J u — ioo i — 1 

(A7) 

where s = u + t' . 

For A'' ^ 1 we can evaluate this integral by the method of stationary phase. Denoting the quantity in square brackets as 
g{s,t), the dominant contribution to the integral comes from near the points (so, to) at which dg/ds = dg/dt = 0, that is, 

JV-l JV-l 

C^T ^' . (A8) 

^ S0+ toAi ^ S0+ toAi 

z — 1 I— 1 

We now show that sq and to are real. Since E, C, u, and A^ are real, the imaginary part of these equations reads 

wo9(so) + wi9(to) = 0, wi^{so) + W2^{to) = 0, (A9) 

where 

N-l .„ 

= r./A ' . D{M,so,to)^mso) + mo)Mf + lQ{so) + Q{to)A,f. (AlO) 
D(Ai,so,to) 

Thus either 9(so) = 9(to) = or 

-^l^Lyy r (All) 

2 fr'^ D{Ai,so,to)D{Aj,so,to) ^ ' 

must vanish. This condition is only satisfied in the trivial case when all of the eigenvalues A; are equal. Thus so and to must 
be real. 

Let the real numbers Sc and tc denote the locations where the integration contours cross the real axis in the complex s 
and t planes. In equation (|A7|I Sc — u and tc = 0. In the method of stationary phase the contours are deformed to cross the 
real axes at s^ = so and tc = to- During this deformation Sc and tc cannot cross the poles at s + tA; = 0. Each such pole 
defines a line in the {sc,tc) plane and in its journey from (u, 0) to (so, to) the point {sc,tc) cannot cross any of these lines. 
This constraint implies that 

t>---^ ifs<0; *>-T^ ifs^O (A12) 

-^^min -'^max 

where Amin is the smallest eigenvalue (other than Ao = 0) and Amax is the largest, i.e., = Ao < Amin ^ Ai, i = 1, . . . , N ^ 

Amax . 

To evaluate the integral (|A7|) we expand the argument g{s,t) of the exponential in a Taylor series around (so,to). We 

have 

^(so,to) = «.o, a^Kto) = ^i, — (so,to) = ^2, where ^„ = ^ (A13) 

i — l 

which is consistent with (lAlOP since so and to are now know to be real. Then 



,iV-3 /■ioo 



cjjv (-E,C)~ — / dt I dsexp 



zoo J u — loo 



g{so,to) + ^wo{s — sq)^ + wi{s — so)(t — to) + iw2(t — to)^ 



(A14) 



We introduce new integration variables {v,v') hy s — sq + iv, t = to + i{v' — vwi/w2) and the integral becomes 

lon{E,C)~ j-^ I duexpl V \ I dv exp[--W2{v) ] 

N-l 

= , = exp toB + soC- Vlog(so+toA,) ; (A15) 

the argument of the square root is positive, as seen from equation ([Alip . 
The entropy of the microcanonical ensemble is 

JV-l ^ 

S{N, E, C) = loga;iv(£', C) = [to-E + sqC - ^ log(so + toAi)j - - log(ij;oW2 - ™?) + const, (A16) 

i— i 

where Boltzmann's constant has been set to unity. For fixed values of so, to. A; the content of the square bracket grows linearly 
with TV while the second term grows only logarithmically. Thus as A'' — > oo we may drop the second term, and 

JV-l 

S{N, E, C) = toE + soC-y log(so + toA,) + const. (A17) 
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The temperature T is defined by 

1 f dS\ f^dto ^dso\ 1 fdso , dto\ 



Using equations (|A8|I . this and the analogous expression for dS/dC simplify to 



-to, (^) =so. (A19) 



We now investigate whether the solution to equations (|A8P exists and is unique. By manipulating the identity A'' — 1 = 
X^ilii^C^o + ioAi)/(so +toAi) it is easy to see that 

soC + toE = A - 1. (A20) 
With this result, equations ([ASP can be combined to give 

r^^^h{N,C,r-to) where h{N, C, T; t) ^ A^A - 1 + tC(A. - F)]-^ ^ 

This non-linear equation can be solved for to given E, C, and A. Furthermore we have 



dt^ ' ' 2[Et7^[A-i + tc(A.-F)]-f ,fe[^-i + ^^(^-r)F[^-i + ^^(^^-r)F' ^ ^ 

which is negative-definite, so h decreases monotonically with t. Thus there is at most one solution for to for given E, C, and 
A. Using equation (|A20|I to eliminate so from the constraints (IA12|) . and recalling that Amin ^ F ^ Amax (see the discussion 
following eq. IA2|I we find that to is restricted to the range 

A- 1 A- 1 

= " C(A^ax-F) < '° < = C(F-A^i„) • 

Note that tmin < and tmax > 0, and 

h{t 

max) — -^min ■ 

(A24) 

Thus h{t) decreases monotonically from Amax to Amin as t grows over its allowed range from tmin to tmax- Since F must lie 
between Amax to Amin for any initial condition, equation (|A2H) always has a unique solution for to given E, C, and A, and 
equation (|A20|) then gives a unique solution for so- 

Note that h(N, C,T;0) = X^fcTi^ A^/A = Aa, the arithmetic mean of the eigenvalues. Thus if F > Aa, to < (negative 
temperature) while if F < Aa then to > (positive temperature). 

The one-particle distribution function fi{Pi,Qi) is defined so that f[P\,Q\)dP\dQ\ is the probability that Pi, Q\ lie in 
the phase-space volume element dPidQi. We have 

„JV-1 N-l JV-1 

fi{Pi,Qi) = — TTT^ / n '^P^'^Q^ ^ \^ - + 0?) - E ^^(P? + s\c^p?-Qi-J2 + 

ojn{E,C)J l\ I ^ J L ^ J 

iJN-1 \e - Ai (P2 + Ql),C- Pi - Qf\ 
= ^ TTrjT. ^- (A25) 

Since A ^ 1 the differences between the numerator and denominator in the solutions of equations (|A8|) for so and to are 
negligible, as are the differences in Wk defined by equation (|A13|) . With this approximation and the use of equation (|A15|) we 
have 

/i(Pi, Qi) = £i±Mi exp [ - (so + toAi)(Pi' + 0?)] . (A26) 

Similar arguments show that the two-particle distribution function is the product of one-particle functions, /2(Pi, Qi, P2, Q2) = 
f{Pi,Qi)f{P2,Q2), etc. 

The average energy in a single normal mode is 

E, = K{Pl + Ql)= A ■ (A27) 

So + toAi 

where the last equality follows from (IA26I) . If so = each mode contains the same energy; in other words, there is equipartition 
of energy or the power spectrum is independent of the frequency Ai. We shall call this 'white noise' - a small abuse of language 
since in most applications white noise corresponds to constant energy per unit frequency rather than per mode. If so < 
then Ei grows as the frequency A^ decreases (red noise). If so > the energy per mode grows as the frequency grows (blue 
noise). The condition so = corresponds to F = A/ X^ilii^ A^^ = A^f , the harmonic mean of the eigenvalues. Since Ah < Aa 
(harmonic mean is less than the arithmetic mean), we have three cases, 

(i) Positive temperature, red noise: 
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(TV - l)A,^in {N -1) TV - 1 

(ii) Positive temperature, blue noise: 

T. A TV- 1 TV- 1 

Aif<r<AA, < So < ^ , — -; — >to>0; 

O OA// 

(iii) Negative temperature, blue noise: 

- < - ^ < » < '^^y <»"• > -c&J- 

Microcanonical ensembles characterized by red noise tend to exhibit warps, because there is more power in low- frequency, 
long-wavelength modes. 



Al Simulations of the microcanonical ensemble 

We want to generate realizations of a relaxed disc containing TV stars on circular orbits, with masses and semi-major axes rm, 
Oi, i = 0, . . . , TV — 1. The disc is specified initially by its energy E and angular-momentum deficit C. The procedure is: (i) 
compute the matrices A (eg. I4U|I and O (eg. I42[) and the eigenvalues A^ of A. (ii) Set F = E/C; if F lies outside the interval 
[Amin, Amax] Spanned by the non-zero eigenvalues Ai, i = 1, TV — 1 then these values of E and C cannot be achieved by a disc 
with the given masses and semi-major axes. Otherwise, find the unique solution to to the non-linear equation (|A21|) and then 
find So from the linear equation (|A20|) . (iii) Set Pq = Qo = 0. Choose (Pi, Qi),j = l,...,TV — fat random from the Gaussian 
distribution (|A26|) . (iv) Set pi = 12^=1 OijPj, Qi = OijQj, then find the inclinations Ii and nodes Hi of the TV stars 

using equations p6[) . The resulting disc will have energy E and angular-momentum deficit C that are within 0(TV~^/'^) of 
the assumed initial conditions (although the differences can still be substantial if the energy or angular-momentum deficit is 
dominated by a small number of modes). 

Note that discs with different values of E and C but the same value of F = E/C have to oc 1/C (by eq. IA2ip and 
So oc 1/C (by eg. IA20|) so the distributions of Pi, Qi, pi, qi, and E scale as -^C (by eg. IA26p . Thus, apart from this trivial 
scaling, the relaxed discs with a given distribution of stellar masses and semi-major axes form a one-parameter family defined 
by the frequency parameter F. 

We generate microcanonical ensembles of discs having the same distribution of stellar masses and semi-major axes as 
in the simulations described at the start of ^ The eigenvalues A; in this disc span a range of 10^ in magnitude (Fig. |3]). 
The largest eigenvalues arise because the interactions of stars that happen to have similar semi- major axes, Aa/a <C 1, 
are unrealistically strong, essentially because the approximation that the Hamiltonian is quadratic is only valid so long as 
the inclinations satisfy / <C Aa/a. Thus the calculations presented here may not accurately represent the thickness of the 
equilibrium disc. 

A sample simulation with frequency parameter F = (ICyr)"^ is shown in Figure \KT\ The two plots show /cos 57 and 
7sinf2 as a function of semi-major axis; the visible high-mass stars (m > 20 M0) are marked with solid blue circles, and the 
low-mass stars are marked by open red circles. The vertical axis is arbitrary since the disc properties can be scaled to other 
energies, as described above. 

As discussed above, the high-frequency modes are unphysical, because the interactions between adjacent stars are un- 
realistically strong. However, most of the structure visible in Figure \KT\ arises from low-frequency modes, which are largely 
independent of the strong interactions between close neighbors. To illustrate this, we have experimented with softening the 
interaction potential according to equation (|37|l . We find that softenings as large as e = 0.003 have very little effect on the 
appearance of realizations of the disc. 

One of the striking features of the disc shown in Figure lAll is the overall warp. To characterize the warp amplitude 
quantitatively, we first define the inclination vectors I; = /i(cos fii, sin f2i) where E and Qi are the inclination and node of 
star i. Then we define the mean and dispersion of the Na inclination vectors in some semi-major axis range A by 

U-^El^. ^3. = 5(^E(I-I.f; (A29) 

the factor of two in the second equation arises because a^. is the dispersion in one of the two components of the vector I. In 
the small-angle approximation within which we are working, the angle between the mean orbit normals in regions A and A' 
is 

e=\lA-lA'\- (A30) 

A measure of the ratio of the warp to the disc thickness is then 6 /a where — ^{ca + i^a')- This ratio is independent of 
the energy E and angular-momentum deficit C so long as the frequency parameter F = E/C is fixed. 
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Figure Al. A realization of the microcanonical ensemble. The upper and lower panels show /cosf! and /sinf2 for a disc with 500 stars 
in the mass range 1 to 30 Mq; the mass function is given by equation I I14I I with a = —0.45 and the semi-major axis distribution is given 
by equation |(T]| with 7 = —1.4. The inner and outer semi-major axes are 1 arcsec and 15.5 arcsec. Solid blue circles denote stars with 
mass > 20 M0, corresponding roughly to stars that are visible in current surveys, and open red circles denote stars with smaller masses. 
The frequency parameter F of eq. I)A1|| equals (10^ yr)~^. 




0.01 1 

10-6 10-5 0.0001 0.001 0.01 0.1 

r(yr->) 

Figure A2. Dimensionless disc warping 9/a as a function of the frequency parameter F. The stellar semi-major axes and masses are 
chosen by the same procedure used to produce Fig. lAll with a different pseudorandom number seed for each point. The parameter 9 is 
the angle between the mean inclination vectors in the inner and outer third of the disc, and c is a measure of the thickness of the disc 
(see text). Magenta and green points correspond to simulations with softening e = and 0.003, respectively. Open circles and crosses 
represent discs with red and blue noise (F less or greater than the harmonic mean of the eigenvalues, Au). Red discs exhibit large warps 
and blue discs do not. 
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Warps are associated with large-scale, low- frequency normal modes, and therefore are stronger in microcanonical ensem- 
bles dominated by red noise. This is illustrated in Figure E2l which shows the warp ratio 6 /a for the visible stars (m > 20 M0) 
in a set of disc realizations with semi-major axes and stellar masses chosen as described at the start of this subsection. Ma- 
genta points are from unsoftened simulations and green points are from simulations with softening e = 0.003. The two regions 
compared, A and A' , are the inner and outer third of the disc stars. Discs with red and blue noise (F < Ah and F > Ah, 
respectively) are marked by open circles and crosses. The warp ratio 6/g declines from ~ 10 for F ^ Ah (substantial warp) 
to ~ 0.3 for F ^ Ah (negligible warp). The softened and unsoftened simulations exhibit the same behavior, except that 
realizations with F>0.003 yr~^ are not present in the softened simulations because the very high-frequency normal modes 
involving two adjacent stars are suppressed by the softening. We conclude that the microcanonical ensemble can exhibit 
significant large-scale warps, depending on the initial conditions as described by F. 
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